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I show that the power iteration method applied to the phase retrieval problem converges under 
special conditions. One is given the relative phases between small non-overlapping groups of pixels 
of a recorded intensity pattern, but no information on the phase between the groups of pixels. 
Numerical tests show that the inverse block iteration recovers the solution in 1 iteration. 



I. INTRODUCTION 

Given a set of intensity measurements, a 2 G M M , an 
unknown object ip G represented by a complex n x 
n image (N = n 2 ), a known "illumination matrix" or 
support matrix Q (C MxAr matrix), a known propagation 
operator (typically one, or a stack of 2D FFT operators) 
F of dimension M x M and set of frames z G C M , which 
are related by: 

z = FQ^, \z\=a. 

Our goal is to find ip or the intermediate variable z, given 
F, Q and a. To do so, we need to find a phase <p such 
that z = a<j) is in the range of FQ. 

We can eliminate ip by using the operator Pq to project 
a vector z onto the range of .FQ: 

P Q = FQ(Q*Q)" 1 Q*F*. (1) 

Which ensures that the unknown vector ip can be ob- 
tained from the frame z by i/j = (Q*Q) _1 Q*F*z. 
A popular approach is to find a vector z such that: 

||[/-P Q ]z|| 2 = 0, (2) 
||[7-P a ]z|| 2 = |||z|-a|| 2 = 0, (3) 

are satisfied simultaneously, and where the Fourier mag- 
nitude projection P a when applied to a vector z, yields: 

P a z = a t^t , a = Diag (a) (4) 
l z l 

where division are intended as element-wise 
operations^'. 



II. PHASE OPTIMIZATION 

Here we want to minimize Eq. [2]w.r.t. a phase vector 
4> = 1, Vi). That is, we want to find: 

argmin \\[I - P Q ]Diag(a)0|| 2 , 

argmin qb*a [I — Pq] a</>, (5) 

<f> 

I discuss three approaches that relax the phase modu- 
lus condition (0*0$ = 1 Vi) to synchronize the relative 
phases. 



a. Power iteration By changing variable z = a</>, we 
write: 

argmin z* (I — Pq) z (6) 

z 

By relaxing (af |^| 2 = a 2 , Vi) and using ||z|| = \\acj)\\ 2 = 
||a|| 2 , we can re- write Eq. (pi as finding the eigenvector 
with largest eigenvalue 3 . Since ||z|| = ||a0|| = ||a|| is 
constant, we rewrite Eq, ([5| as: 

argmax z*Pqz, (7) 

z 

we apply one step of power iteration: 

v t+1 = Pqz (8) 

We then form a projection on the unit torus to ensure 
that = 1 (or \zi\ = ai) by element- wise normaliza- 

tion: 

M 

Here we have obtained the classical alternating projection 
method, which is known to stagnate with classical CDI 
but to converge (slowly) in ptychographic imaging. 

b. Greedy phase optimization Since the diagonal 
term z*Diag(PQ)z is also independent on the choice of 
<j) (for \(j)i\ = 1), one can remove it when computing the 
power iteration: 

= [P Q -Diag(P )]zW (9) 

After we apply the projection of to the unit torus, 
we obtain the following updated 

=P a (P Q - Diag (Pq))z^ 

In classical CDI, Pq.. = j^p- is simply the sum of the 
support volume (or area) normalized by the oversampled 
volume, in ptychographic imaging Diag(PQ) is the ratio 

of intensities Pq.. = |^np for every pixel of a frame i 
generate by a submatrix Qi- At the first iteration, using 
data generated from the object in Fig. [TJwith a random 
phase as a starting guess, Eq. Q appears to out-perform 
Eq. (|6|, however the two methods converge to similar lo- 
cal minimum within ten iterations. The relaxations in 
Eqs. ( |6|9[ ) are similar. By removing diagonal compo- 
nents we change the relaxation. In Eq. ^ we have 
||a 2 0|| constant, in Eq. ^ ||[1 - Diag (Pq)]ct</>|| is con- 
stant. However Diag (Pq) is often constant and the two 
relaxations are equivalent, giving more weight to high 
intensity values. 
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c. 



Inverse iteration 1 . If we solve the minimization 
problem (Eq. ([5| with a different relaxation, setting 
1 1 <j)\ | = k to a constant, we re- write the problem as 



arg max 



H 



i[I-P]t 



(10) 



and apply the power iteration: 



This method is commonly referred to as inverse iteration 
and it is used to find the smallest eigenvector of a matrix. 
We note however that any v written in the following way: 



(a _1 P Q a) x 



(12) 



is an eigenvector with eigenvalue of a(7 — Pq)a, there- 
fore the inverse iteration method cannot be applied di- 
rectly. 

When H is singular, then instead of power iteration we 
may want to find the smallest modification of the phase 
that is in the null space of which we can write it as a 
re- weighted LSQ problem of the form: 



arg mm 



(z - J) 



S.t. Z = PqZ 



(13) 



which provides a search direction toward the solution 
that differs from standard projection algorithms. An- 
other approach is to include additional restrictions on v 
before applying the inverse iteration as described in the 
following. 

d. Inverse block iteration it was observed that 



computing the exact solution to Eq. ( 11 ) after "binning' 



or fixing the relative phase between groups of pixels, 
improved convergence rate in large scale ptychographic 
imaging. The use eigensolvers for the interferometric case 
was also suggested iiP, for the connection Laplacian of a 
graph. 

Let us introduce a binning matrix T* composed of a 
series of masks T{ that integrate over a region of dimen- 
sion M/k of the data (in Fourier domain). For example, 
we can partition our data in 3, creating a tall matrix of 
dimension M x 3: 



1m/3> Om/3 5 Om/3 



T = I M / 3 , l M / 3 , M / 3 
Om/3: Om/3? l-M/3 



where T\ 



V A M/3' u M/3' 



0^y 3 )* is a vector of length M. 



We restrict our search of the solution to Eq. (11) by 
restricting v to be: 



(14) 
in Eq. 

(11) we obtain the inverse iteration step with initial 0- 



v = Diag ((jf) Too, 
If we multiply from the left by T*Diag I 



phase vector as first guess: 



hV+ 1 ) = a 1 i, 



(15) 



Where 



= T*Diag (>>)* [/ - P Q ]Diag (z^) T, 

Ai is a scalar multiplicative factor, and 1 is a vector of 
appropriate length (= 3 in this example). By comput- 
ing co>^ +1 ) from Eq. (15), and from Eq. (14), and 



(11) projecting on the unit torus we obtain the update 



Diag 



\Tco( e + 1 )\ 



(16) 



In the following section we'll show an example of the 
inverse iteration method. 



III. NUMERICAL EXAMPLE 

Here ip consists of the cameraman image of 32 x 32 
pixels, embedded in a matrix of 64 x 64 pixels (Fig. [I]). 
The "illumination matrix" is the support of the object, 
Q = Diag (S). The support is 1 inside the 32 x 32 box 
containing the image, and otherwise. (Q*Q) _1 is re- 
placed by the pseudoinverse Q = Q*. The Fourier trans- 
form of if; was perturbed by 32 x 32 randomly distributed 
phases (Fig. |3|, each multiplying a bin of 2 x 2 pixels 
(Fig. J2|. Upon perturbation, the image in real space 
(Fig. HT cannot be distinguished. Many iterations of Eq. 
[6] or Eq. [9] cannot converge (Fig. [5] show ing Eq. ([6] ) 
updates), while 1 iteration of Eqs. ( 15|16 ) converges to 
the solution (Fig. pi). 



IV. CONCLUSIONS 

I have shown that power iteration methods can recover 
phase perturbations under special circumstances. If one 
is given the relative phases between a small group of pix- 
els (binned) and a random perturbation of the phase be- 
tween all the groups of pixels (the bins), then the in- 
verse block iteration can recover the solution in 1 iter- 
ation, it was observed that the inverse block iter- 
ation improved convergence rate in large scale ptycho- 
graphic imaging. The inverse block iteration was also 
shown to recover perturbations in the experimental ge- 
ometry such as position errors and intensity fluctuations. 
More work is needed to determine the optimal combina- 
tion of Eqs. ( 6|9|13|15|16 ), and the properties of T, in 
large scale phase retrieval problems. 
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DISCLAIMERS 
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sponsored by the United States Government. While 
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or implied, or assumes any legal responsibility for the ac- 
curacy, completeness, or usefulness of any information, 
apparatus, product, or process disclosed, or represents 



that its use would not infringe privately owned rights. 
Reference herein to any specific commercial product, pro- 
cess, or service by its trade name, trademark, manufac- 
turer, or otherwise, does not necessarily constitute or im- 
ply its endorsement, recommendation, or favoring by the 
United States Government or any agency thereof, or the 
Regents of the University of California. The views and 
opinions of authors expressed herein do not necessarily 
state or reflect those of the United States Government or 
any agency thereof or the Regents of the University of 
California. 
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FIG. 1. Object ip (64 x 64) used to simulate diffraction 
data 
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FIG. 3. Random perturbation: (64 x 64) phases gener- 
ated by 32 x 32 random phases each spread over a bin of 
(2 x 2) pixels. 
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FIG. 2. Each column of the matrix T extracts an area 
of 2 x 2 pixels out of an image of (64 x 64) pixels. 




FIG. 4. Image in real space (|.F*z|) after random phase 
perturbation (Fig. [3| 





FIG. 5. Image in real space after random phase pertur- 
bation (Fig. |j), using Eq. [^updates ( \F* (P a P s ) 1000 z\). 



FIG. 6. Image in real space after one step of Eqs. ( 15|16 ) 
update. 



